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ABSTRACT 

We present numerical simulations of the growth and saturation of the Kelvin- 
Helmholtz instability in a compressible fluid layer with and without a weak mag- 
netic field. In the absence of a magnetic field, the instability generates a single 
eddy which flattens the velocity profile, stabilizing it against further perturba- 
tions. Adding a weak magnetic field - weak in the sense that it has almost no 
effect on the linear instability - leads to a complex flow morphology driven by 
MHD forces and to enhanced broadening of the layer, due to Maxwell stresses. 
We corroborate earlier studies which showed that magnetic fields destroy the 
large scale eddy structure through periodic cycles of windup and resistive decay, 
but we show that the rate of decay decreases with decreasing plasma resistivity 
7], at least within the range of t] accessible to our simulations. Magnetization 
increases the efficiency of momentum transport, and the transport increases with 
decreasing rj. 

Subject headings: instabilities — methods: numerical — MHD — turbulence 

1. INTRODUCTION 

Turbulent mixing layers, found at the boundary of two fluids in relative motion, are 
thought to entrain ambient material into outflows and to contribute to the near uniform 
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chemically mixed interstellar medium (ISM) (see reviews by Scalo & Elmegreen 2004 and 
Elmegreen & Scalo 2004). Shear flow instabilities (ie. Kelvin Helmholtz Instability (KHI)) 
provide the main mechanism of forming these mixing layers. Because shear flows show 
up in many places in astrophysics and geophysics, the properties of KHI has been studied 
extensively using both linear stability analysis and non linear simulations (see §2.2p . 

In this paper, we expand on these results, focusing primarily on the effects of a weak 
magnetic fleld, too weak to affect the initial development of the instability. A scenario for 
the nonlinear evolution of weakly magnetized shear layers was proposed by Frank et al. 1996 
and Malagoli et al. 1996, based on numerical simulations. The fleld is wound up in the flow 
until the point of reconnection, which injects energy into smaller scale eddies in the flow. The 
fleld then winds up in these eddies until it reconnects again. This process is continued until 
the perturbed magnetic energy is damped out, leaving an enlarged shear layer. It remains to 
be shown that this picture holds for resistivities t] appropriate to the ISM, which are much 
lower than anything that can be achieved in a numerical simulation. In this paper, we study 
the evolution of the layer with t] varied by a factor of 50. We flnd that certain properties - 
notably the time for the instability to saturate, the energy at saturation, and the energy put 
into the magnetic fleld - appear to saturate with r]. Others - the rate at which perturbations 
decay and the rate at which the layer broadens with time - do not reach r] independent states 
for the range of rj accessible in our calculations. 

Section [2] describes the basic problem setup and discusses the linear theory. Details on 
the numerical method and parameters of our models are given in ^ The results (§!]) are 
summarized in ^ 



The conservative variables that describe the system are the density p, momentum density 
pv, total energy density pE, and magnetic induction B. The evolution of these variables is 
given by the MHD equations. Written in divergence form, these are: 
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f + V . (vB - Bv) = 'J^V'B, (4) 

where we have included terms for both viscous and resistive dissipation. In the case of an 
incompressible fluid and a spatially isotropic viscosity, the viscous term in the momentum 
equation reduces to V-II = p/xV^v, where fi is the coefficient of viscosity. Numerical schemes 
cannot resolve the dissipation scales of the ISM, which is very close to the ideal limit, with 
a Reynolds number. Re = ^ ~ 10^'' and a magnetic Reynolds number, = ^ ~ 10^^. 
We will discuss our treatment of viscous and resistive dissipation in § 13. 2[ The gas pressure 
p is related to the energy density by p = (7 — l){pE — |pv^ — ^) where 7 is the adiabatic 
gas constant, which we keep fixed at 5/3. 

The equilibrium fiow is aligned along the x direction and is sheared in the y direction, 

vo = x^tanh(-), (5) 
2 a 

while the gas pressure and density are initially constant, p = po, p = Pq. Eq. [5] describes a 
fiow that is in the positive x direction for y > and in the negative x direction for y < 
and smoothly varies within a shear layer of total width 2a, asymptotically reaching ±^ for 
III ^ 1. In the absence of perturbations, the velocity profile will broaden due to viscous 
dissipation on a timescale t^sc- In practice, we restrict ourselves to timescales t such that 
t -C tyisc- In the MHD case, we take the equilibrium field to be spatially constant and aligned 
with the fiow: 

Bo = xSo. (6) 

It will be convenient to describe the problem with dimensionless parameters. The ratio 
of the fiow speed and the sound speed, Cs^ = is the sound Mach number, Mg = 

The ratio of the fiow speed to the Alfven speed, Cafl = is the Alfven Mach number. 

Ma = The strength of the magnetic field can be described by the ratio of the Alfven 
speed to the sound speed, a = ^ = The level of viscous diffusion is characterized by 
the Reynolds number. Re = — and that of resistive diffusion by the magnetic Reynolds 
number Rm = where L is the characteristic length of the system. 



2.2. Linear Stability 

A linear stability analysis will tell us under what conditions an instability will arise 
and how fast we expect it to develop. A linear stability analysis involves perturbing the 
linearised MHD equations and solving for the growthrate F of the instability as a function 
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of the perturbed wavenumber, k. The hnear properties of shear flow instabihties have been 
extensively studied elsewhere. What follows is a summary of salient results. 

Chandrasekhar (1961) solved the incompressible, inviscid, ideal vortex sheet problem 
where the velocity profile is characterized by a discontinuity: Vq = for y > 0, Vq — 
for y < 0. In the purely hydrodynamic case (Bq = 0), assuming p is the same in the two 
media, the growthrate is Thd = |k • When the flow velocity and wave vector line up, 
all modes are unstable with the smallest scales growing the fastest. 

A magnetic fleld that is oriented perpendicularly to the flow plane has no effect on the 
stability of the flow. It only acts to make the fluid more incompressible. However, a magnetic 
field aligned parallel to the fiow has a stabilizing effect. The tension in the magnetic field 
tends to suppress the instability, modifying the growthrate to: 



If the initial velocity, magnetic field, and wave vector all line up, the condition for stability 
is Ma = ^ < 2. In other words, the Alfven speed has to be larger than half the shear 
velocity difference (ca,o > -y) to stabilize the layer. 

The effects of a shear layer, from a linear profile (eg. Ray 1982) to a tanh profile (Lau 
& Li 1980, Miura & Pritchctt 1982) have also been looked at. Ferrari & Trussoni (1983) 
summarized the previous work by examining profiles that range between these two cases. 
The effects of different profiles are all qualitatively the same, though the exact details depend 
on the initial equilibrium. In general, introducing a layer adds another lengthscale to the 
problem. For wavelengths that are long compared to the size of the layer, the details of the 
layer become unimportant and the growthrate approaches the vortex sheet result. However, 
when the wavelength is short compared to the layer width, 1/k ^ a, the perturbation no 
longer "sees" the shear and becomes stable. As a result, the fastest growing modes are no 
longer the shortest wavelengths but are now ones whose wavelength is comparable to the 
size of the layer, ka ~ 1. A parallel magnetic field still tends to stabilize the fiow, though 
the stability criterion (M^ > 2) only holds approximately. The exact criterion depends on 
the details of the problem. 

The effect of compressibility on shear fiows has been studied by Landau (1944), Miles 
(1957), Fejer (1963) and others for a vortex sheet, and by Blumen (1970), Blumen, Drazin, 
& Billings (1975), Drazin & Davey (1977), Ray (1982), and Ferrari & Trussoni (1983) for 
a shear layer. Compressibihty has the general effect of lowering the growthrate from the 
incompressible case. Also, Miles (1957) showed that compressibihty adds an upper Mach 
number limit on the stability for a vortex sheet for parallel perturbations (it remains unstable 
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to 3D perturbations): for Ms > V2 the flow becomes stable. However, in the presence of a 
shear layer, all Mach numbers are unstable, but the growthrates in the supersonic case are 
about an order of magnitude less than in the subsonic case. The spatial distribution of mode 
energy is sensitive to Ms as well: for subsonic flow the modes decay exponentially away from 
the layer, while in supersonic flow they decay much more slowly (Blumen et al. 1975). 

Hughes & Tobias (2001) summarize and expand the literature regarding the differences 
between a 2D and 3D treatment of the linear theory of shear flow instabilities. In a purely 
hydrodynamic flow, for every unstable mode in 3D, there is a corresponding unstable mode in 
2D at a lower Reynolds number. Re = aVo/fi. Therefore, as Re is increased to a value that is 
no longer stabilizing, the 2D mode will become unstable first (Squires 1933) As a consequence, 
a 2D stability analysis of hydrodynamic flows is adequate. Michael (1953)and Stuart (1954) 
showed that, for an MHD flow, for every unstable 3D mode there is a corresponding 2D mode 
at both a lower Reynolds number Re and magnetic Reynolds number R^- In this case, the 
region of neutral stability is a function of both Re and Rm- Therefore, the 2D mode will 
not necessarily become unstable first (Hunt 1966). In the limit of ideal MHD {Re = oo and 
Rm = oo) an unstable 2D mode is a faster growing mode than the corresponding 3D mode 
(Hughes & Tobias 2001). Therefore, although it is adequate to study ideal MHD flows in 
2D, non-ideal MHD flows may have to be treated in 3D. 

In order to develop some feeling for the effects of resistivity and viscosity, and for the 
effect of a 3D wavevector for the set of models we ran, we carried out a parameter study 
using a linear eigenvalue code. Our code solves the linearized incompressible MHD equations 
with respect to a background flow and magnetic field. The linear perturbation is Fourier 
decomposed in the x and z directions, and the remaining y direction is discretized with 
a Chebyshev spectral method (Boyd 2001, Trefethen 2000). The eigenvalue problem is 
solved in an infinite domain y G [—oo, oo], which is mapped into y' G [—1, 1] by a mapping 
scheme that localizes most of the grid points near y = 0. This ensures that we resolve the 
instability, and is compatible with the boundary conditions that all perturbations vanish as 
y goes to ±oo. We used 150 grid points along y to perform a parameter scan in viscosity 
and resistivity to get a handle on the relative effects on both. We have also looked at the 
effect of an increasing component to the perturbation. In both studies, we have used the 
equilibrium flow described in § 12.11 and initial parameters given in § 13. 1[ 

Figure [1] shows contours of growthrate as a function of both the viscous, (/i), and 
resistive, {i]), diffusivity. As long as the viscous diffusivity is ;^ 0.0003, the growth rate is 
insensitive to resistivity and decreases with increasing viscosity, eventually becoming stable. 
At low viscous diffusivity (yU ^ 0.0003) the effect of viscosity becomes negligible. Resistivity 
only starts to make a non-negligible difference for rj ^ 0.001, above which it acts as a 
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destabilizing mechanism. These two observations make sense because viscosity acts as a 
dissipation mechanism for the kinetic energy of the instabihty and resistivity, though it is a 
dissipation mechanism for the magnetic energy, also affects the frozen-in level of the magnetic 
field, which is a stabilizing mechanism. As rj increases the field decouples from the fluid and 
the growthrate approaches the hydrodynamic growthrate. As t] decreases, the field becomes 
more frozen-in and will then approach the ideal MHD growthrate. Since we are considering 
a weak field case only, the amount of stabilization provided by the magnetic field is small 
and therefore the differences between the MHD and hydro growthrate is small. 

Fig. [2] shows the growthrate as a function of kz with a constant kx for similar flow 
properties as given in Table [H The fastest growing mode occurs when kz = and decreases 
for growing k^ until it becomes stable. Since we are considering a weak-field problem in which 
resistivity and viscosity have little effect on the growthrate (see Fig{T]), it is not surprising 
that the 2D mode is the most unstable. 



3. THE COMPUTATIONS 
3.1. Models Run 

We focus on the non-linear evolution of the Kelvin-Helmholtz instability in weak-field 
MHD flow. For realistic values of viscosity and resistivity, the instability growthrates in 
the two flows are similar (see Fig. [T]). However, the nonlinear evolution of the two are very 
different (Frank et al. 1996, Malagoli et al. 1996). As explained in the Introduction, we wish 
to examine the role of resistivity in the evolved flow, and to compare momentum transport 
in the HD and weak fleld MHD flows. 

Table [T] provides a summary of the models run. The numbers in the names of the 
MHD runs refer to the magnetic Reynolds number Rm- Even though the maximum value 
Rm = 50000 is far below that of the ISM [R m ~ 10 ) we are able to draw useful conclusions 
(see §4.3). The entries in the last column denote "low" (256x512), "medium" (512x1024), 
"high" (1024x2048), and "super high" (2048x4096) resolution, respectively; the multiple 
runs being used to test convergence. All the runs have the same initial pressure (po = 1-0), 
density (po = 1-0), sonic Mach number Mg = 1 and the MHD runs all have Ma = 10, so 
a = 0.1 and (3 = 120. We solve eqns. (1) - (4), with the equilibrium flow given by eqn. [5] 
and eqn. [6l and an initial velocity perturbation given by: 



6Yy = 6Vo sm{kxx) exp(— 




where the amplitude SVq <^ The attenuation with distance from the layer on a lengthscale 
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given by a ensures that the perturbation stays locahzed to the y = interface. 

The computations are done on a 2D grid in the x-y plane. The dimensions of the grid 
are x = [0, L] and y = [—L, L] where L is the characteristic length of the system. We 
choose the characteristic length such that a single wavelength fits in the box, L = 27r/A;^.. 
We have chosen kx = 27t such that the growthrate is near a maximum (see Keppens et al. 
2001), and we set a/L = 0.05 and a/L = 0.2. These are the same initial conditions used 
by both Malagoli et al. (1996) and Keppens et al. 2001). These choices allow a small, yet 
resolvable shear layer, and an attenuation scale that is both larger than the layer (a/a = 4) 
and shorter than the characteristic length, so that the ±y boundaries interact minimally 
with the instability. We have chosen the characteristic time to be the sound crossing time 
tg = L/csfi. All the models presented here are run over a timescale of 20t<j, which is long 
enough to capture the saturation and at least the initial non-linear evolution of the instability. 

3.2. Numerical Method 

The magnetohydrodynamical scheme - Proteus - is based on a conservative gas-kinetic 
fiux-splitting method, introduced by Xu (1999) and Tang & Xu (2000). Viscosity and ohmic 
resistivity (RHS of eq. S]) are implemented via dissipative fluxes (Heitsch et al. 2004, 2007). 
Thus, Proteus allows the full control of dissipative effects. Test cases have been presented 
by Heitsch et al. (2007). Proteus integrates equations ([Hll]) at 2nd order in time and space. 
We use only the resistivity implementation for the models presented here, relying on the 
(extremely low) numerical viscosity. 

Proteus is also equipped with a Lagrangian tracer mechanism (Heitsch et al. 2006), 
which permits us to follow particles advected by the flow and thus enables us to study 
mixing. We defer a discussion of the mixing properties of the flow to a future paper (Palotti 
et al. 2007) except for a brief mention in §4.1. 

4. RESULTS 
4.1. Convergence 

The numerical dissipation scale is related to the resolution (although usually not di- 
rectly analytically to Ax^/At). As resolution increases, the dissipation scale decreases, and 
structure is allowed to develop on those smaller scales. 

Setting the dissipative scales allows us to reach detailed convergence (i.e. a "point- 
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to-point" agreement in physical variables between resolutions) in the linear regime of the 
instability. Detailed convergence usually cannot be reached in the non-linear (or turbulent) 
regime, however, virtual convergence (i.e. the convergence in integrated quantities between 
resolutions, such as magnetic energy) is still possible if the problem is nearly independent of 
the dissipation scales. We will be using "convergence" in the latter, virtual, sense. 

In a weak field MHD flow, the magnetic field becomes wound up in the flow, increasing 
the current density and enhancing the role of resistive dissipation. As we show in the ensuing 
discussion, resistive relaxation is rather abrupi0. The resistive event alters the structure of 
the flow. Therefore, resistivity plays an integral part in determining the evolution of the 
system. Viscosity, on the other hand, becomes important when energy cascades down to the 
dissipation scale and appears to be unimportant in determining the bulk properties of the 
flow. For these reasons, we use a physical resistivity while relying on numerical viscosity. 

We have carried out a convergence study by running each of our models (see Table [1] at 
resolutions ranging from 256x512 to 2048x4096. In Fig. [3] we present the (virtual) conver- 
gence of the total kinetic energy in the y direction, KEy = J J ^pVydxdy, as a function of time 

for the HD model, as well as the total perturbed magnetic energy, AME = J J — ^:^dxdy, 
for the models MHD5000 and MHD50000. 

Any discrepancies in the different resolutions for the KEy in the HD model can be 
attributed to numerical viscosity only. It appears that KEy is a converged quantity. Even 
though there is not true convergence, the differences between the 512 x 1024 and 1024 x 
2048 models are less than 5% while those of the 256x512 and 512x 1024 can be as high as 
10%, ie. the differences are small and decrease with increasing resolution. 

The perturbed magnetic energy is a good indicator of our ability to resolve the resistive 
length scales. For the larger resistivity, = 5000, is well converged. For times t < IGtg 
the difference between the 512x1024 and 1024x2048 runs is less than 2% and reaches at 
most ~ 5.5%. Achieving convergence in the MHD50000 model is much more difficult. For 
times t < lAtg the differences between the 1024x2048 and 2048x4096 runs are less than ~ 
5%. At later times, the differences reach up to ~ 15%. This is an indication that we have 
not fully resolved the resistive length scale. However, the rate of dissipation of the magnetic 
energy qualitatively follows the same trend for both the 1024x2048 and 2048x4096 runs. 
Therefore, for t < lAtg, which includes both the saturation of kinetic energy (t = GAtg) 
and the saturation of the magnetic energy {t = 9.0ts), all the runs are effectively converged 



^We do not have the resolution to see whether the flow resembles any of the standard modes of recon- 
nection, such as Sweet-Parker reconncction, but it has been shown elsewhere that shear flow instability can 
enhance the reconnection rate (Knoll & Chacon 2002). 
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and for later times, runs with Rm <^ 5000 are effectively converged while those with lower 
resistivity only show qualitatively similar trends. 

As discussed in §2.1, we chose the initial horizontal velocity to be proportional to 
tanh (y/a) (see eqn. ([5])), with the parameter a, which we call the layer width, set to 0.05. If 
we take horizontal averages of Vx at constant y at later times, (see Fig. [10] for examples), the 
transition region between Vq and —Vq is well fit by a straight line. From this we can define 
the width of the layer as the position where the velocity, v = We find the layer width is 
an increasing function of time. In Figure H] we present the velocity layer width as a function 
of time for the HD run and the MHD run MHD5000. Over the entire time run, both layer 
profiles agree to within a few percent for the highest resolution runs. This, together with 
the convergence of KEy and AME, is good evidence that the global properties of the flow 
are well converged. 

As we mentioned in §3.2, we can trace the trajectories of tracer particles (for the higher 
resolution runs, however, this requires a prohibitively large amount of computer time). In 
order to test the convergence of the small scale structure, we computed the separation of 
pairs of particles as a function of time for the run MHD 1000. The two highest resolution 
runs agree for times less than 15ts. The differences after this can be attributed to the effects 
of numerical viscosity on the flow. We defer other discussion of the tracer particles to Palotti 
et al. (2007). 

For timescales t <^ IQts, each of the models can be considered converged. The lack of 
convergence at longer times is a result of the development of small scale flow structure and 
consequent viscous dissipation, which is entirely numerical. As the resistivity increases, less 
small scale structure develops, and thus the convergence is better. At the lowest resistivities, 
numerical resistivity may be becoming important as well. 

4.2. Nonlinear Flow Morphology 

The growth rates in the HD and MHD models differ very little, because the initial 
magnetic field is very weak (see Table 1). Over the range of resistivities we examined 
{Rm = 1000 — 50000), we calculated the growth rates to differ by ~1%, increasing as the 
resistivity increases. The trend reflects the progressively weaker coupling of the field to the 
fluid with increasing resistivity, and hence a reduction of the stabilizing magnetic tension 
force. 

However, the nonlinear evolution of these flows are vastly different. In the HD flow, a 
single large eddy develops and remains until viscosity eventually damps it out, on timescales 



much longer than considered here. The total kinetic energy in the y-direction, the top panel 
in Fig. [3], is consistent with a tumbling elliptical eddy. The amplitude of each successive 
oscillation decreases due to viscosity. The density structure and velocity field at time t = SAtg 
for the HD model are plotted in Fig. O The density shows a rarefied central region around 
which the single eddy spins. The density contrast between the central region and the outer 
region is ~ 36% in this case. The tumbling eddy is also evident in the velocity field (see 



The evolution of the MHD model is much more complicated, and certain aspects of 
it depend on resistivity. We have plotted the magnetic energy structure (grey scale) and 
magnetic field (vectors) for three times, 6.4ts,12ts, and 20^^, and two values of Rm, 5000 and 
50000, in Figure These times represent the saturation of KEy {GAtg), the end of the 
run (20ts) where, at Rm ^ 5000, both KEy and AME have nearly decayed fully, and an 
intermediate time (12ts) where the energies have partially decayed. Our primary observations 
about Fig. are that (1) the magnetic field is wound up by the large scale eddy, and that 
this wound up field induces a wealth of magnetic structure, that (2) reducing the resistivity 
allows the structure to exist at smaller scales, and that (3) reducing the resistivity allows 
the structure to persist for longer times. Some of the magnetic folds and loops seem to be 
magnetic islands, which can only form through resistive evolution. Whether this is better 
described as reconnection or diffusion is unclear, however (see footnote 1). Along with the 
magnetic structure, there is a corresponding small-scale density and velocity structure in the 
MHD runs, with density perturbations about ± 2% for R^, = 5000 and a maximum of about 
2% and minimum of about 9% for R^, = 50000. 



Saturation of the KHI is defined as the point at which the kinetic energy in the y 
direction peaks in magnitude. We found the time of kinetic energy saturation in both the 
HD and all MHD runs to be 6.4^^ (see top panel of Fig. [3] and Fig. [T^ . The mechanism for 
saturation is the interplay between the y components of the forces on the plasma, namely gas 
pressure, Reynolds stress and Maxwell stress. Integrating the y component of the momentum 
equation (eql2]) from x = to x = L and from y = to y = L (the upper half of the domain) 
and applying periodicity in x, we get 



where we have assumed the viscosity to be negligible. When d{pVy)/dt = the y-forces 
balance and the instability is saturated. We have plotted the the various contributors to the 



FiglSD. 



4.3. Saturation and Effect of Resistivity 
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y component of the the force for our MHD50000 and HD runs in FigjTl The Reynolds stress 
term {—pVy) is the component of force that drives the instabihty while the gas pressure {—p) 
component acts to slow down the instability in the HD model. In the MHD models, the 
Maxwell stress term {-^) also acts to slow down the instability initially, but it is clear from 
Fig 17] that magnetic forces are much smaller in magnitude than hydrodynamic forces in this 
flow. 

We also examined how resistivity affects the saturation. In Fig. [S] we plot the kinetic 
energy in the y direction, the total perturbed magnetic field energy, and the sum of the 
two at tx as a function of magnetic Reynolds number. For the sake of comparison, we 
have included the HD results at Rm = 0, assuming that as the resistivity becomes infinite, 
the MHD flow approaches the HD limit. At t^, the kinetic energy is greatest for the HD 
model and decreases for increasing Rm- The magnetic energy, on the other hand increases 
for increasing magnetic Reynolds number so that the total perturbed energy increases with 
Rm- By the time of saturation, the magnetic energy has increased between 35% (lowest 
Rm) to 48% (highest Rm) over the background value. As Rm increases {Rm > 5000), the 
perturbed kinetic and magnetic energies seem to approach equipartition, with a ratio of 
KEy/AMEtot =~ 1.14, whereas at Rm < 5000, the ratio is 2.8. 

After the kinetic energy saturates, the magnetic energy continues to grow (see Fig. [TTj) . 
eventually saturating at a later time, Im- Both tM and the value of the magnetic energy at 
tM depend on resistivity. We plot both of these in Fig. O As the magnetic Reynolds number 
increases, the time at which the magnetic energy saturates plateaus at tm ~ 9.0, while the 
maximum magnetic energy continues to increase. 

In order to get a better handle on the mechanism behind the magnetic saturation, we 
calculated the volume integrated time rate of change of the magnetic energy, ignoring the 
surface terms, 

-<^-'''^«»' 

The first term on the right hand side represents the ohmic dissipation while the second term 
represents the work done by field on the gas. We found that for all models run, during the 
saturation and decline of the magnetic energy, the work term is positive meaning that the 
fiow is doing work on the field. Therefore, the decline in magnetic energy is not dynamical - 
the field never becomes strong enough to unwind itself - but is purely dissipative in nature, 
at least in a global sense, over the range in Rm considered here. The increases in tM and 
AME with Rm reflect the fact that as time progresses, the flow puts energy into the fleld 
at smaller and smaller scales. Roughly speaking, the ampliflcation is quenched when the 
energy reaches the resistive scale; this happens sooner, and at lower magnetic energy, for the 
smaller Rm- However, it seems likely that AME is bounded above by the kinetic energy in 
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the flow. The flattening of AME with Rm seen in Fig. ([8]) is evidence for this. Since this is 
a dynamical bound, not a resistive one, we speculate that there is a transition from resistive 
to dynamical saturation at resistivities lower than what we can examine. 

4.4. Time Evolution 

Figure ffTTl) depicts a cycle of growth and decay in both kinetic and magnetic energy for 
Rm = 5000 (similar evolution is observed for smaller Rm as well). After the kinetic energy 
saturates, the magnetic energy continues to grow. Since the growth is at the expense of 
kinetic energy, KEy decreases during this time. As the field grows its structure becomes more 
complex. Eventually the rate of work done by the fiow is overcome by resistive dissipation. 
During the decay of the field the kinetic energy grows again, energized in part by the small 
scale forces characteristic of magnetic reconnection and in part by the background velocity 
shear that drives instability in the first place. The cycle repeats itself, but at lower amplitude 
due to resistive dissipation. After the second peaks, the perturbations simply decay. This is 
the behavior found by Frank et al. (1996) and Malagoli et al. (1996). 

Figure f|T2|) shows that the picture is more complex, and depends on Rm- At Rm = 20000 
and 50000, the decay is much less pronouncecj^. Whereas for Rm < 5000 the perturbations 
have dissipated almost completely by the end of the run, for Rm = 50000 the perturbed 
kinetic and magnetic energy densities are still more than half their maximum values. The 
temporal evolution is less cyclic, and shows bumps and dips that might be associated with 
discrete magnetic reconnection or dynamical relaxation events, superimposed on a general 
decay. 

Most importantly, there is no sign that the decay rate has reached an asymptotic limit 
independent of Rm- It seems to us that there are two possibilities for the large Rm behavior. 
One is that the decay rate is determined by a dissipation coefficient other than resistivity. 
Since the viscous diffusivity exceeds the magnetic diffusivity in most of the ISM, viscosity 
is a likely candidate. Another possibility, if the gas is weakly ionized, is ion-neutral friction. 
We are currently simulating fiows in which ion-neutral damping dominates; Palotti et al. 
(2008). The alternative is that the decay rate is determined by a turbulent diffusivity that 
is independent of any microscopic diffusion coefficient and is attained at higher Rm than we 



^ Since the computations at R,n = 50000 arc not as well converged as for the lower R,n (see Fig. [3]), we 
have plotted values of the energy for the highest resolved run in this case as asterisks. Comparison of the 
asterisks with the curves only underscores the basic point as the decrease in magnetic energy is, if anything, 
slower at the higher resolution 
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can achieve. We have no basis for determining which of these scenarios would hold in the 
ISM, where Rm is ~ 11 orders of magnitude larger than it is in our simulations. 

One way or another, the energy in the layer must eventually be dissipated. Does it 
matter how? We argue that it does, both for the observational appearance of the layer and 
for the transport taking place within it. This is demonstrated in the following subsection. 



4.5. Momentum Transport 

As the instability evolves and KEy increases, momentum is transported between the 
upper {y > 0) and lower {y < 0) domains. As a result, the initial sheared velocity profile 
(see eq. [5]) will spread out. In Fig. [10] we plot the x-averaged x velocity as a function of y at 
different times for both the HD and MHD5000 models. 

The evolution of the profile differs between the two models. In the HD model, the profile 
reaches a maximum width of ~ 0.12 at saturation, t = GAtg (see dashed line in Fig. [TOj) . and 
then the width oscillates in a non-uniform fashion (see Fig. H]), with a minimum of ~ 0.08 at 
(t = 12ts). The average width is ~ 0.1. If we take the evolved profile as an initial condition, 
we find that it is stable (presumably it is unstable to longer wavelength modes that are not 
captured within our periodic domain). 

The profile width in the MHD model, however, continues to grow after saturation, finally 
reaching its maximum of 0.37 at t ~ ISt^, when the perturbations have nearly dissipated, 
and remaining constant thereafter. The spatial roughness of the MHD case, which contrasts 
with the smooth profile in the HD case, is due to the prominent small scale flow structure in 
the MHD case. Note that only a portion of the domain is plotted in each figure; the shear 
layer has not reached the top and bottom boundaries of the domain, and the momentum 
flux escaping through these boundaries is found to be small. 

In order to understand why the two models differ, we average the x component of the 
momentum equation (eql2]) over x, getting 

^ = ^{{Pv.vy) - (B^By)) (11) 

where the first term on the RHS represents Reynolds stress and the second term represents 
Maxwell stress. We plot each of these components for the MHD model at various times 
in Fig. [131 Initially, the Reynolds term dominates the spreading of the layer, reaching a 
maximum value at saturation [t = 6.4^^) and then starts to decrease in magnitude, following 
the HD model. However, the Maxwell stress term takes over and further broadens the layer 
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until t ~ 15 when the broadening ceases. In both cases, the speed at which the layer broadens 
is subsonic. 

Because in the MHD case the layer is broadened by turbulent stresses, and the amplitude 
and decay rate of the turbulence depend on i?^, we expect the evolution of layer width to 
depend on Rm as well. This is borne out by Fig. (fT4|) . The layer width increases with Rm, 
and, like the decay rate, shows no sign of leveling off as Rm increases. Thus, momentum is 
transported more efficiently in an MHD model compared to an HD model, and the efficiency 
increases with increasing Rm- Our range of Rm is too small and too far from ISM values to 
reliably extrapolate. 

5. SUMMARY AND CONCLUSIONS 

We have run 2D simulations exploring the non-linear evolution of the Kelvin Helmholtz 
Instability in both a magnetized and unmagnetized shear flow. We concentrated on two 
models: a hydrodynamic model at sound Mach number, Ms = 1, and a weak- field magneto- 
hydrodynamic model with sound and Alfven Mach numbers, Mg = 1, = 10. We have con- 
sidered the effect of resistivity by running models with Rm ranging from Rm = 1000 — 50000, 
at resolutions from 256 x 512 to 2048 x 4096. 

Justification for restricting the problem to two dimensions was provided in part by our 
parameter study of the linear regime: using a linear eigenmode code, we showed that the 
2D growthrate is the fastest growing mode for the parameters we considered. However, the 
growthrate does not give any information as to the non-linear evolution of the flow, so it is 
possible that 3D effects would become important in the nonlinear regime. 

The linear growth phase and the saturation mechanism of a HD flow and a weak field 
MHD fiow show many similarities. We calculated the growthrates to be different by only 
~ 2.5% between the HD and MHD50000 models. The mechanism for saturation is an 
interplay between the Reynolds stress trying to spread the layer and the combined effort of 
gas pressure and Maxwell stress (in MHD only) trying to prevent the spreading (see Fig. [7]). 
We found that the gas pressure dominates over the Maxwell stress in preventing the spreading 
of the layer. However, as the initial magnetic field increases, the contribution of Maxwell 
stress will become more important in saturating the instability, and, for a sufficiently large 
magnetic field, will suppress it entirely. 

Despite these similarities in the initial phase of growth, the non-linear evolution between 
the two models is drastically different (see §4.2p . as has been known for some time (Frank 
et al. 1996, Malagoli et al. 1996). The differences are evident in both the evolution of 
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the kinetic and magnetic energies (top panel of Fig. [3] for HD and Fig. [TT] for MHD) 
and in the HD density structure (Figs. [5]) and MHD magnetic structure (Fig. [6]). The 
HD flow develops a single large tumbling eddy that remains for the duration of the runs 
(viscosity will eventually damp it out). In the MHD flow, the eddy winds up the magnetic 
field. Magnetic forces spawn a rich array of small scale flow, density, and magnetic field 
structures. As the resistivity is lowered, the energy in the fluctuations increases (although 
it is probably bounded by the equipartition value) and the decay time of the fluctuations 
increases. Previous works have suggested that resistive dissipation is the dominant process 
in the evolution of the flow. However, if the decay time of the fluctions continue to increase 
for resistivites lower than considered here, a process with a larger diffusion coefficient may 
become the dominant process in controlling the evolution of the energy in the flow. In 
the ISM, processes with larger diffusion coefficients include both viscosity and ion-neutral 
friction. Including these processes is beyond the scope of this work. 

The decay rate of the turbulence, and the mode of decay, is not just of academic interest. 
We showed that the shear layer is broadened by turbulent Maxwell stresses well beyond what 
occurs in the hydrodynamic case. The longer the turbulence survives, and the larger its 
amplitude, the broader the layer becomes. This means that momentum transport increases 
with decreasing resistivity. 

Although the turbulent decay rate and the momentum transport rate are sensitive to 
Rm over the range accessible to our study, other properties of the instability are relatively 
insensitive. The instability growth rate, times of saturation, and saturated amplitude are 
relatively independent of Rm, primarily because the magnetic field is too weak to affect these 
properties. Nevertheless, the dependence on Rm found for some of the important quantities 
suggest that extrapolation of these results to real astrophysical flows should be selective and 
requires some caution. 

The authors are happy to acknowledge support from NSF grants AST0507367 and 
PHY0215581 to the University of Wisconsin and NASA grant NNG06GJ32G to the Univer- 
sity of Michigan. Portions of this paper were written at the KITP at UC Santa Barbara; we 
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Table 1. Models Run 



Model 




Ma 


a 






Res 


HD 


1 


OO 





oo 





l,m,h 


MHDIOOO 


1 


10 


0.1 


120 


1000 


l,m 


MHD2000 


1 


10 


0.1 


120 


2000 


l,m 


MHD5000 


1 


10 


0.1 


120 


5000 


l,m,li 


MHD20000 


1 


10 


0.1 


120 


20000 


l,m,h,sli 


MHD50000 


1 


10 


0.1 


120 


50000 


l,m,h,sh 



- 19 - 




-6 -S -4 -3 -2 

log 0(1 



Fig. 1. — Contours of growthrate are plotted as a function of the viscous, rj, and resistive, 
/i, diffusivities for a flow speed and magnetic field comparable to those given in Table [TJ 
Viscosity acts to stabilize the instability by dissipating kinetic energy while the resistivity 
acts to destabilize the instability by decoupling the field from the flow. At low viscosities, 
the growthrate only varies by a few percent because we are considering only a weak magnetic 
field. 
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Fig. 2. — The growthrate as a function of is plotted for flow speeds and magnetic fields 
comparable to those given in Table [H The fastest growing mode for the parameters con- 
sidered is the kz = mode. As a consequence a 2D study should be an adequate place to 
start. 
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Fig. 3. — The total kinetic energy in the y direction as a function of time for the HD model 
(top) and the perturbed magnetic energy as functions of time for model MHD5000 (middle) 
and MHD50000 (bottom) are plotted for three resolutions. 




Fig. 4. — The width of the velocity profile as a function of time for both the HD (top) and 
MHD5000 (bottom) models. 
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Fig. 5. — Gas density for the high resolution hydro dynamic model after saturation. The 
velocity vector field is overplotted. 
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Fig. 6. — Logarithm of the magnetic energy for the high resolution Rm— 5000 (left) and 
i?^=50000 (right) at three times after saturation. The magnetic vector field is overplotted. 
The times shown are 6.4, 12, and 20 tg. 




Fig. 7. — The total Reynolds stress, gas pressure, and Maxwell stress as functions of time 
for the HD (top) and MHD50000 1024x2048 resolution (bottom) models. Saturation occurs 
when all the forces cancel each other out. 
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Fig. 8. — The total y kinetic energy, perturbed magnetic energy and total perturbed energy 
at saturation as a function of resistivity. 
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Fig. 9. — The maximum pcrtm^bcd magnetic energy (top) and the time to reach the maxi- 
mum perturbed magnetic energy (bottom) as a function of resistivity. 
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Fig. 10. — The x-averaged x-velocity as a function of y at various times for botii tiie HD 
(top) and MHD5000 (bottom) models. 
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Fig. 12. — The y-componcnt of the kinetic energy, KEy (top) and the perturbed magnetic 
energy AME (bottom) as a function of time for magnetic Reynolds number ranging from 
Rjn — 1000 to Rjn — 50000. Because the — 50000 run is not as well converged, we have 
also included the highest resolution data as the asterisks. 



-30- 




Fig. 13. — The Reynolds and Maxwell stress as a function of y for six times at Rm = 5000. 
The times represented are 6 (top) and 12 (bottom) tg. 
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Fig. 14. — The width of the velocity profile as a function of time for magnetic Reynolds 
number ranging from Rm — 1000 to Rm — 50000. 



